home *** CD-ROM | disk | FTP | other *** search
/ Atari Mega Archive 1 / Atari Mega Archive - Volume 1.iso / gnu / gawk / gawk213s.zoo / gawk-src-2.13 / atof.c < prev    next >
C/C++ Source or Header  |  1991-05-27  |  13KB  |  570 lines
  1. /*
  2.  *     double strtod (str, endptr);
  3.  *     const char *str;
  4.  *     char **endptr;
  5.  *        if !NULL, on return, points to char in str where conv. stopped
  6.  *
  7.  *     double atof (str)
  8.  *     const char *str;
  9.  *
  10.  * recognizes:
  11.  [spaces] [sign] digits [ [.] [ [digits] [ [e|E|d|D] [space|sign] [int][F|L]]]]
  12.  *
  13.  * returns:
  14.  *    the number
  15.  *        on overflow: HUGE_VAL and errno = ERANGE
  16.  *        on underflow: -HUGE_VAL and errno = ERANGE
  17.  *
  18.  * bugs:
  19.  *   naive over/underflow detection
  20.  *
  21.  *    ++jrb bammi@dsrgsun.ces.cwru.edu
  22.  *
  23.  * 07/29/89:
  24.  *    ok, before you beat me over the head, here is a hopefully
  25.  *    much improved one, with proper regard for rounding/precision etc
  26.  *    (had to read up on everything i did'nt want to know in K. V2)
  27.  *    the old naive coding is at the end bracketed by ifdef __OLD__.
  28.  *    modeled after peter housels posting on comp.os.minix.
  29.  *    thanks peter!
  30.  *        ++jrb
  31.  */
  32. #if !(defined(unix) || defined(minix))
  33. #include <stddef.h>
  34. #include <stdlib.h>
  35. #include <float.h>
  36. #endif
  37. #include <errno.h>
  38. #include <assert.h>
  39. #include <math.h>
  40.  
  41. #ifdef minix
  42. #include "lib.h"
  43. #endif
  44.  
  45. #ifndef _COMPILER_H
  46. #include <compiler.h>
  47. #endif
  48.  
  49. extern int errno;
  50. #if (defined(unix) || defined(minix))
  51. #define NULL     ((void *)0)
  52. #endif
  53.  
  54. #define Ise(c)        ((c == 'e') || (c == 'E') || (c == 'd') || (c == 'D'))
  55. #define Isdigit(c)    ((c <= '9') && (c >= '0'))
  56. #define Isspace(c)    ((c == ' ') || (c == '\t'))
  57. #define Issign(c)    ((c == '-') || (c == '+'))
  58. #define IsValidTrail(c) ((c == 'F') || (c == 'L'))
  59. #define Val(c)        ((c - '0'))
  60.  
  61. #define MAXDOUBLE    DBL_MAX
  62. #define MINDOUBLE    DBL_MIN
  63.  
  64. #define MAXF  1.797693134862316
  65. #define MINF  2.225073858507201
  66. #define MAXE  308
  67. #define MINE  (-308)
  68.  
  69. /* another alias for ieee double */
  70. struct ldouble {
  71.     unsigned long hi, lo;
  72. };
  73.  
  74. static int __ten_mul __PROTO((double *acc, int digit));
  75. static double __adjust __PROTO((double *acc, int dexp, int sign));
  76. static double __ten_pow __PROTO((double r, int e));
  77.  
  78. /*
  79.  * mul 64 bit accumulator by 10 and add digit
  80.  * algorithm:
  81.  *    10x = 2( 4x + x ) == ( x<<2 + x) << 1
  82.  *    result = 10x + digit
  83.  */
  84. static int __ten_mul(acc, digit)
  85. double *acc;
  86. int digit;
  87. {
  88.     register unsigned long d0, d1, d2, d3;
  89.     register          short i;
  90.     
  91.     d2 = d0 = ((struct ldouble *)acc)->hi;
  92.     d3 = d1 = ((struct ldouble *)acc)->lo;
  93.  
  94.     /* check possibility of overflow */
  95. /*    if( (d0 & 0x0FFF0000L) >= 0x0ccc0000L ) */
  96. /*    if( (d0 & 0x70000000L) != 0 ) */
  97.     if( (d0 & 0xF0000000L) != 0 )
  98.     /* report overflow somehow */
  99.     return 1;
  100.     
  101.     /* 10acc == 2(4acc + acc) */
  102.     for(i = 0; i < 2; i++)
  103.     {  /* 4acc == ((acc) << 2) */
  104.     asm volatile("    lsll    #1,%1;
  105.              roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  106.         : "=d" (d0), "=d" (d1)
  107.         : "0"  (d0), "1"  (d1) );
  108.     }
  109.  
  110.     /* 4acc + acc */
  111.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  112.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  113.  
  114.     /* (4acc + acc) << 1 */
  115.     asm volatile("  lsll    #1,%1;
  116.              roxll   #1,%0"    /* shift L 64 bit acc 1bit */
  117.         : "=d" (d0), "=d" (d1)
  118.         : "0"  (d0), "1"  (d1) );
  119.  
  120.     /* add in digit */
  121.     d2 = 0;
  122.     d3 = digit;
  123.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  124.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  125.  
  126.  
  127.     /* stuff result back into acc */
  128.     ((struct ldouble *)acc)->hi = d0;
  129.     ((struct ldouble *)acc)->lo = d1;
  130.  
  131.     return 0;    /* no overflow */
  132. }
  133.  
  134. #include "flonum.h"
  135.  
  136. static double __adjust(acc, dexp, sign)
  137. double *acc;    /* the 64 bit accumulator */
  138. int     dexp;    /* decimal exponent       */
  139. int    sign;    /* sign flag          */
  140. {
  141.     register unsigned long  d0, d1, d2, d3;
  142.     register          short i;
  143.     register           short bexp = 0; /* binary exponent */
  144.     unsigned short tmp[4];
  145.     double r;
  146.  
  147. #ifdef __STDC__
  148.     double __normdf( double d, int exp, int sign, int rbits );
  149.     double ldexp(double d, int exp);
  150. #else
  151.     extern double __normdf();
  152.     extern double ldexp();
  153. #endif    
  154.     d0 = ((struct ldouble *)acc)->hi;
  155.     d1 = ((struct ldouble *)acc)->lo;
  156.     while(dexp != 0)
  157.     {    /* something to do */
  158.     if(dexp > 0)
  159.     { /* need to scale up by mul */
  160.         while(d0 > 429496729 ) /* 2**31 / 5 */
  161.         {    /* possibility of overflow, div/2 */
  162.         asm volatile(" lsrl    #1,%1;
  163.                     roxrl    #1,%0"
  164.             : "=d" (d1), "=d" (d0)
  165.             : "0"  (d1), "1"  (d0));
  166.         bexp++;
  167.         }
  168.         /* acc * 10 == 2(4acc + acc) */
  169.         d2 = d0;
  170.         d3 = d1;
  171.         for(i = 0; i < 2; i++)
  172.         {  /* 4acc == ((acc) << 2) */
  173.         asm volatile("    lsll    #1,%1;
  174.                  roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  175.                  : "=d" (d0), "=d" (d1)
  176.                  : "0"  (d0), "1"  (d1) );
  177.         }
  178.  
  179.         /* 4acc + acc */
  180.         asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  181.         asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  182.  
  183.         /* (4acc + acc) << 1 */
  184.         bexp++;    /* increment exponent to effectively acc * 10 */
  185.         dexp--;
  186.     }
  187.     else /* (dexp < 0) */
  188.     { /* scale down by 10 */
  189.         while((d0 & 0xC0000000L) == 0)
  190.         {    /* preserve prec by ensuring upper bits are set before div */
  191.         asm volatile(" lsll  #1,%1;
  192.                     roxll #1,%0" /* shift L to move bits up */
  193.             : "=d" (d0), "=d" (d1)
  194.             : "0"  (d0), "1"  (d1) );
  195.         bexp--;    /* compensate for the shift */
  196.         }
  197.         /* acc/10 = acc/5/2 */
  198.         *((unsigned long *)&tmp[0]) = d0;
  199.         *((unsigned long *)&tmp[2]) = d1;
  200.         d2 = (unsigned long)tmp[0];
  201.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  202.         tmp[0] = (unsigned short)d2;    /* the quotient only */
  203.         for(i = 1; i < 4; i++)
  204.         {
  205.         d2 = (d2 & 0xFFFF0000L) | (unsigned long)tmp[i]; /* REM|next */
  206.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  207.         tmp[i] = (unsigned short)d2;
  208.         }
  209.         d0 = *((unsigned long *)&tmp[0]);
  210.         d1 = *((unsigned long *)&tmp[2]);
  211.         /* acc/2 */
  212.         bexp--;    /* div/2 taken care of by decrementing binary exp */
  213.         dexp++;
  214.     }
  215.     }
  216.     
  217.     /* stuff the result back into acc */
  218.     ((struct ldouble *)acc)->hi = d0;
  219.     ((struct ldouble *)acc)->lo = d1;
  220.  
  221.     /* normalize it */
  222.     r = __normdf( *acc, ((0x03ff - 1) + 53), (sign)? -1 : 0, 0 );
  223.     /* now shove in the binary exponent */
  224.     return ldexp(r, bexp);
  225. }
  226.  
  227. /* flags */
  228. #define SIGN    0x01
  229. #define ESIGN    0x02
  230. #define DECP    0x04
  231. #define CONVF    0x08
  232.  
  233. double strtod (s, endptr)
  234. const register char *s;
  235. char **endptr;
  236. {
  237.     double         accum = 0.0;
  238.     register short flags = 0;
  239.     register short texp  = 0;
  240.     register short e     = 0;
  241. /*    const register char *start = s;  */
  242.     double       zero = 0.0;
  243.  
  244.     assert ((s != NULL));
  245.     
  246.     if(endptr != NULL) *endptr = (char *)s;
  247.     while(Isspace(*s)) s++;
  248.     if(*s == '\0')
  249.     {    /* just leading spaces */
  250.     return zero;
  251.     }
  252.     
  253.     
  254.     if(Issign(*s))
  255.     {
  256.     if(*s == '-') flags = SIGN;
  257.     if(*++s == '\0')
  258.     {   /* "+|-" : should be an error ? */
  259.         return zero;
  260.     }
  261.     }
  262.     
  263.     do {
  264.     if (Isdigit(*s))
  265.     {
  266.         flags |= CONVF;
  267.         if( __ten_mul(&accum, Val(*s)) ) texp++;
  268.         if(flags & DECP) texp--;
  269.     }
  270.     else if(*s == '.')
  271.     {
  272.         if (flags & DECP)  /* second decimal point */
  273.         break;
  274.         flags |= DECP;
  275.     }
  276.     else
  277.         break;
  278.     s++;
  279.     } while (1);
  280.  
  281.     if(Ise(*s))
  282.     {
  283.     if(*++s != '\0') /* skip e|E|d|D */
  284.     {  /* ! ([s]xxx[.[yyy]]e)  */
  285.     
  286.         while(Isspace(*s)) s++; /* Ansi allows spaces after e */
  287.         if(*s != '\0')
  288.         { /*  ! ([s]xxx[.[yyy]]e[space])  */
  289.     
  290.         if(Issign(*s))
  291.             if(*s++ == '-') flags |= ESIGN;
  292.     
  293.         if(*s != '\0')
  294.         { /*  ! ([s]xxx[.[yyy]]e[s])  -- error?? */
  295.  
  296.             for(; Isdigit(*s); s++)
  297.             e = (((e << 2) + e) << 1) + Val(*s);
  298.  
  299.             if(IsValidTrail(*s)) s++;
  300.             
  301.             /* dont care what comes after this */
  302.             if(flags & ESIGN)
  303.             texp -= e;
  304.             else
  305.             texp += e;
  306.         }
  307.         }
  308.     }
  309.     }
  310.  
  311.     if((endptr != NULL) && (flags & CONVF))
  312.     *endptr = (char *) s;
  313.     if(accum == zero) return zero;
  314.     
  315.     return __adjust(&accum, (int)texp, (int)(flags & SIGN));
  316. }
  317.  
  318. double atof(s)
  319. const char *s;
  320. {
  321. #ifdef __OLD__
  322.     extern double strtod();
  323. #endif
  324.     return strtod(s, (char **)NULL);
  325. }
  326.  
  327. #ifdef TEST
  328. #ifdef __MSHORT__
  329. #error "please run this test in 32 bit int mode"
  330. #endif
  331.  
  332. #define NTEST 10000L
  333.  
  334. #ifdef unix
  335. #ifdef __MSHORT__
  336. #define    RAND_MAX    (0x7FFF)    /* maximum value from rand() */
  337. #else
  338. #define    RAND_MAX    (0x7